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We study a matrix product state (MPS) algorithm to approximate excited states of translationally 
invariant quantum spin systems with periodic boundary conditions. By means of a momentum 
eigenstate ansatz generalizing the one of Ostlund and Rommer pQ, we separate the Hilbert space 
of the system into subspaces with different momentum. This gives rise to a direct sum of effective 
Hamiltonians, each one corresponding to a different momentum, and we determine their spectrum by 
solving a generalized eigenvalue equation. Surprisingly, many branches of the dispersion relation are 
approximated to a very good precision. We benchmark the accuracy of the algorithm by comparison 
with the exact solutions of the quantum Ising and the antiferromagnetic Heisenberg spin-1/2 model. 
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I. INTRODUCTION 

Recently we have presented an algorithm [2] for the 
approximation of the ground state of translationally in- 
variant (TI) spin chains with periodic boundary condi- 
tions (PBC) my means of TI MPS. In this work we will 
use the ground states obtained in [2] as the basis of an 
ansatz for excited states with definite momentum. Wc 
will consider only spin chain Hamiltonians that are trans- 
lationally invariant thereby fulfilling [H, T] = where 
T is the translation operator that shifts the lattice by 
one site. Furthermore, as we will deal with finite chains 
in the following, it means that there is no well defined 
momentum operator for our systems. Nevertheless we 
can classify translationally invariant states by their quasi- 
momentum which is defined in terms of the their eigen- 
value with respect to T. This definition is sensible since 
in the thermodynamic limit, if wc keep the chain length 
fixed, the lattice spacing becomes infinitesimally small 
and the quasi-momentum becomes identical to the mo- 
mentum, which is then well defined. For convenience, 
we will use the term momentum when we actually re- 
fer to the quasi-momentum. This should not cause any 
confusion since we will only deal with quasi-momenta 
throughout this work. 

Since H and T commute, they can be diagonalized si- 
multaneously. This suggests that any variational ansatz 
based on eigenstates of the translation operator will be 
well suited to define families of states within which min- 
imization with respect to some variational parameters 
will yield momentum eigenstates with minimal energy. 
Formulating this observation in terms of an MPS-based 
ansatz has led in the past to some very interesting re- 
sults about excitation spectra. The first approach in this 
direction has been made in [1 where the main result is 
the celebrated insight that the fixed point of the density 
matrix renormalization group (DMRG) can be written as 
an MPS. In addition to this, based on the MPS that is 
obtained for the ground state of the infinite Heisenberg 
spin-1 chain, the authors suggest a variational ansatz for 



excitations with definite momentum. Since the trans- 
lationally invariant MPS they start with is an approxi- 
mation of the ground state in the thermodynamic limit, 
their ansatz for excitations is only well suited in the limit 
N — > oo. For finite chains, the idea of using momentum 
eigenstates for the diagonalization of TI Hamlitonians 
has been used in [3] in order to obtain a few of the lowest 
branches of excitations for the bilinear-biquadratic (BB) 
spin-1 chain. The resulting state is a TI superposition 
of a special class of tensor network states, which can be 
viewed as an extension of MPS with PBC [4] to states 
that can accommodate multipartite entanglement. Even 
though the multipartite entanglement is a nice feature 
which yields a better variational ansatz in the cases when 
the approximated states have that special entanglement 
structure (in [3] one has in addition to the usual maxi- 
mally entangled virtual bonds between nearest neighbors 
a virtual GHZ state connecting all sites) we will not adopt 
it in our present ansatz. Furthermore we would like to 
point out that the individual MPS tensors produced by 
the minimization procedure in [3] are not TI, only their 
superposition is. 

Recent results [2] on the approximation of ground 
states of TI PBC Hamiltonians opened up the possibility 
of unifying the ideas from pQ and [5J in order to obtain 
an algorithm for excitations with definite momentum in 
which only one local tensor has to be determined, thereby 
avoiding the usual sweeping procedure and the associated 
factor N in the computational cost. One of the main fea- 
tures of TI MPS is the fact that the tensor network that 
has to be contracted for the computation of expectation 
values contains big powers of a so-called transfer ma- 
trix [5]. For non-critical systems the eigenvalues of this 
transfer matrix usually decay rapidly enough s.t. big 
powers thereof can be accurately approximated by con- 
sidering only a few dominant eigenvectors. In these cases 
the computational cost for the evaluation of expectation 
values for systems with PBC can be reduced significantly 
from 0(D 5 ) to 0(D 3 ), where D denotes the virtual bond 
dimension of the MPS. For critical systems however the 
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eigenvalues of the transfer matrix decay much slower and 
the algorithm that must be employed in order to obtain 
the optimal approximation within the class of MPS with 
fixed D has a scaling that depends in a not yet fully un- 
derstood way [2] on D, N and on the universality class 
of the simulated model. 

The ansatz we present in this work is based on TI- 
MPS and thereby all computed quantities will contain 
big powers of the transfer matrix. We would like to em- 
phasize that the computational cost can be reduced by a 
factor of D 2 only in the case of non-critical systems. For 
critical systems the full contraction of tensor networks 
(i.e. without using any approximations of the transfer 
matrix) will turn out to have a more favorable overall 
scaling of the computational cost. Details on why this 
is the case and on the scaling of the computational cost 
can be found in the next section. 



tually very similar. An important feature of all these ap- 
proaches is the reduction of the dimension of the problem 
by a factor N. This is reached by effectively projecting 
the original problem into the subspace with fixed momen- 
tum k and minimizing the energy within the variational 
manifold spanned by the free parameters in the ansatz. 
In our case these free parameters are the components B 
of an MPS tensor. As it is always the case with MPS 
algorithms, one must eliminate the ambiguities arising 
from the MPS representation by fixing the gauge. Here 
this is done by starting with certain tensors A in ^ 
and not changing them throughout the entire minimiza- 
tion procedure. This automatically fixes the gauge of the 
tensors B as they are surrounded on both sides by A. 

III. THE ALGORITHM 



II. OVERVIEW 

Due to T N — 1, the translation operator T that shifts a 
state on a PBC lattice with N sites by 1 site is the gener- 
ator of the cyclic group of order N. Hence its eigenvalues 
Tk must be roots of the unity i.e. = e~ lh Tr with inte- 
ger fee [0,7V — 1]. An ansatz for eigenstates of T with 
eigenvalue e~ is obviously given by 
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l^(B)) = J2 
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L T"|0 A (B)) 
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Henceforth we will refer to states of the form ([!]) as Bloch 
states. Note that we have used the convention that T is 
the operator that realizes a translation by one site to the 
right s.t. T\4>(i 1 ,i 2 ,...,iN)} = \4>{iN,i\, «2> • • • , ijv-i))- 
The state |</>a(B)) can in principle be any arbitrary state, 
but in order to exploit the advantages of TI MPS, we 
choose 
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with identical matrices on all sites except the first one. 
We will choose the Ai to be the matrices corresponding to 
the best TI MPS ground state approximation for a given 
model. We emphasize that the Ai remain fixed through- 
out the entire simulation. This is the reason why we have 
omitted them from our labeling convention for the Bloch 
states |V'fc(B)). We have used bold letters in order to 
denote objects that are obtained if one rearranges the 
components of three indexed MPS tensors into vectors, 
i.e. A := vec(A i ° /3 ). After fixing the momentum k, the 
Bloch states |V>fc(B)) will depend only on the tensors B 
which will define the variational manifold. 

Our ansatz for Bloch eigenstates differs slightly from 
the ones presented in Refs. [TJ [21 [5] although it is concep- 



Ansatz ([I]) defines a class of variational states for the 
lowest energy states with fixed momentum. The energy is 
a quadratic expression in the tensor B and thereby, as it 
is usually the case in MPS based algorithms, minimizing 

(^.(B)|ff|^(B)) 
& (k) - mm , [6) 

BeC d D 2 (^ fe (B)|^ fe (B)) 

is equivalent to solving a generalized eigenvalue equation 

tf e// (fc)B,(fc) = EiWNtffWBiW (4) 
where H e ff(k) is defined by 



Btff e// (fc)B := (^(B)|Jf|^ fe (B)) 



and N eff (k) by 



BtJV e// (fc)B := (Vfc(B)|^ k (B)) 



(5) 



(6) 



The eigenvector corresponding to the smallest eigen- 
value Eo(k) yields then the tensor Bo (A;) that when 
plugged into our ansatz ([I]) gives the momentum-k state 
with minimal energy. Note that the variational princi- 
ple guarantees that only the Bloch state ([T]) with lowest 
energy is the best approximation to the exact cigenstate 
with that momentum within the subspace spanned by 
our ansatz states. However if the lowest energy state is 
approximated accurately, due to the fact that the other 
Bj(fc) are orthogonal to Bo(fc), the next solution Bi(fe) 
has a good chance to be close to the next higher energy 
state with that momentum. In fact it will turn out that 
quite a few of the higher energy solutions of Q are good 
approximations to low energy states with fixed momen- 
tum. Their precision is most of the time surprisingly good 
given the fact that the variational principle does not hold 
for these states. The quality of these solutions depends 
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FIG. 1. (Color online). Definition of No m (A) as the norm of 
a TI MPS determined by the tensor A. 



strongly on the bond dimension D, the chain length N 
and the model under consideration. 

The bottleneck of our method is the computation of 
the effective matrices H e ff(k) and N e ff(k). Let us first 
consider N e ff(k) since it is the slightly simpler one. ft 
reads 



BtjV e// (A:)B = 
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where 7Vo m (A) is a tensor network resembling the norm 
of a TI MPS with empty slots and m (see figure [IJ . 
To get from the second to the third line we have used 
the fact that due to the PBC only the relative distance 
between n and m plays a role. In the last line we have 
merely renamed the summation index and introduced the 
quantity N 0m (A). Thus in order to obtain N e ff(k) we 
have to compute the contraction of the N tensor net- 
works iVo m (A) and then take the sum of these terms 
after weighting each one of them with the corresponding 
phase factor. The computational cost for the contraction 
of each tensor network is 0(D 6 ) s.t. the overall cost for 
computing N eff (k) is 0(ND 6 ). 



,(A) := 




First, due to the translational invariance of the Hamilto- 
nian we can write 
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(8) 
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where hoi is the term acting between the first two sites of 
the chain. Note that in ^ we have restricted ourselves 
to nearest neighbor Hamiltonians since this is the type of 
models we will treat numerically in this work. General- 
izing the ideas developed here to any local Hamiltonian 
is straightforward. With ph, H e ff(k) reads 
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Again, due to the fact that the rh and n sums run over all 
N sites of a PBC chain, it is irrelevant where they begin 
s.t. the I sum merely yields a factor N. We rename the 
summation indices for convenience and obtain 
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where H 0nm (A) is a tensor network resembling the ex- 
pectation value of an operator acting on the sites n and 
n + 1 with respect to a TI MPS where the slots and 
m have been left open (see figure [2]) . The computational 
cost for the contraction of each tensor network is again 
0(D 6 ) but now we have a total of N 2 summands s.t. the 
overall cost for computing H e ff(k) is 0(N 2 D 6 ). Note 
that to obtain H e tf(jk) is computationally the most ex- 
pensive part of our algorithm so we can say that the 
overall computational cost scales like 0(N 2 D 6 ). 



FIG. 2. (Color online). Definition of Ho nm (A) as the expec- 
tation value of a two-site operator with respect to a TI MPS A - Overall scaling of the computational cost 
determined by the tensor A. 

At first sight the cost seems horrible for a 1D- 
H e ff(k) is constructed very much in the same spirit. algorithm. Let us however have a closer look at what 
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FIG. 3. (Color online). Lowest ten branches of the excitation spectrum for a critical Ising chain with N = 50. Left: D = 8. 
Right: D = 32. 



we get for this price. First of all note that if we com- 
pute the sets of matrices {N 0m (A)} and {H 0nm (A)} for 
n,m € [0, N — 1] and store these, we can obtain the 
H e ff(k) and N e ff{k) for all k trivially by just build- 
ing the appropriately weighted sums. For each of these 
k we then have to solve the generalized eigenvalue equa- 
tion Q. Since H e fAk) and N e ff(k) are small dD 2 x dD 2 
matrices solving Q docs not represent any difficulty 
and can be done using any standard library eigenvalue 
solver. Each eigenvalue problem leads to dD 2 orthonor- 
mal vectors Bj(fc) which plugged into the ansatz (fTl) 
yield dD 2 states. Thus computing the sets {N 0m (A)\ 
and {Ho nm (A)} only once supplies immediately us with 
NdD 2 states! By comparing our numerical results to 
exactly solvable models we will show that the low en- 
ergy states obtained in this way are very accurate. This 
means that in terms of computational time per state our 
algorithm performs quite well. 

The computational bottleneck at the moment is that 
we have to store N 2 dD 2 x dD 2 matrices in the mem- 
ory. With the present implementation, for a chain with 
N = 100 sites, we can go up to D = 32. For larger N 
simulations we have to settle for smaller D. It is how- 
ever straightforward how this boundary can be pushed 
considerably towards larger D. First, instead of keeping 
all matrices in the memory, one can write them to the 
hard disk after computing each of them. Second, since 
the {AT 0m (A)} and {H 0nm (A)} are independent, one can 
parallelize their computation. 

Thus the conceptual bottleneck becomes the contrac- 
tion of the tensor networks {No m (A)} and {Ho nm (A)}. 
For non-critical systems big powers of the transfer matrix 
can be well approximated by a few of its dominant eigen- 
vectors [5] and the contraction of most of the {N 0m (A)} 
can be done with the computational cost 0(n 2 D 3 ) while 
that of most of the {Ho nm (A)} with the cost 0{n 3 D 3 ). 
Here n represents the dimension of the subspace within 
which we approximate the powers of the transfer ma- 
trix [5] . This cannot be done however for critical systems 



where in principle n may grow as big as D 2 thereby yield- 
ing a much worse scaling than the naive 0(D 6 ). Note 
that since {N 0m (A)} and {Ho nm (A)} are open tensor 
networks the 0(D 5 ) contraction scheme [I] that works 
for expectation values (i.e. closed tensor networks) can- 
not be applied here. Additionally, even if we restrict our- 
selves to non-critical systems, not all of the {iVo m (A)} 
and {i?onm(A)} can be computed with the cost that 
scales like D 3 : if the distance between the open slots 
is not big enough, we cannot use the approximation for 
big powers of the transfer matrix between the slots, and 
we are back to exact contraction for this portion of the 
chain which in the case of {iVo m (A)} leads to the over- 
all scaling 0(nD 5 ) and in the case of {H 0nm (A)} to the 
scaling 0(n 2 D 5 ). Thus the very naive exact contraction 
procedure that we use is not so bad after all in this case 
even if it scales like 0(D 6 ). 

There is one more subtlety we would like to point out 
here. It turns out that the matrix N e ff(k) is always 
singular which presents a problem when we try to solve 
the generalized eigenvalue equation ^ since the solu- 
tion involves the inverse N~^(k). We can circumvent 
this problem by solving Q within the nonsingular sub- 
space like it has been done in [T] . Eigenvectors associated 
to the zero eigenspace of N e ff(k) will result in physical 
states IV'fe(B)) = 0, i.e. these are states of zero length in 
the Hilbcrt space. Any physical operator will produce a 
zero when acting on these states. In particular, the ef- 
fective Hamiltonian H e ff(k) will also have zero eigenval- 
ues for the same eigenvectors, and we do not loose any 
information by restricting to the nonsingular subspace. 
The dimension of the zero eigenspace can be shown to 
be D 2 (d - 1) for k ^ and D 2 {d - 1) + 1 for k = as 
we demonstrate in [TJ. The tricky point is that for some 
models the strictly non-zero eigenvalues of N e ff(k) be- 
come so small that they yield the generalized eigenvalue 
problem ill conditioned. In general this behavior does 
not occur for small D. For big D or in certain regions of 
the phase diagram however the nonsingular eigenvalues 
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become so small that it is hard to distinguish numerically 
between the singular subspace and the nonsingular one. 
This issue might be the source of the mysterious nega- 
tive gap that appears in [T] in the vicinity of the critical 
point. 

We have employed a slightly different method for the 
regularization of N e f f(k). Instead of projecting the prob- 
lem into the strictly non-singular subspace, we restrict 
ourselves to the subspace in which the eigenvalues of 
N e ff(k) are larger than some e. There is a tradeoff be- 
tween loss in precision due to this projection and loss in 
precision due to the bad conditioned generalized eigen- 
value problem. In the end we have settled for a seemingly 
optimal e — 10 -11 . 



IV. NUMERICAL RESULTS 

We have applied the algorithm presented above to two 
exactly solvable nearest neighbour interaction spin mod- 
els in order to benchmark its accuracy: the quantum 
Ising model and the Heisenberg spin-1/2 antiferromag- 
netic chain. Even though the Heisenberg model is ex- 
actly solvable by means of Bethe ansatz, in practice it 
is much harder to obtain its entire low-energy spectrum. 
This is due to the fact that the elementary excitations are 
two-spinon states and among these, the solution of the 
Bethe ansatz equations for the two-spinon singlet states 
are computationally very challenging [S]. Thus for long 
chains we have restricted ourselves to check only the pre- 
cision of the lowest two-spinon triplets. For small chains 
that are accessible via exact diagonalization on the other 
hand, we compare not only the entire low-energy spec- 
trum but also the fidelity of the states themselves. 



A. Quantum Ising model 

The Hamiltonian we have used in our simulations of 
the Quantum Ising model is given by 

i i 

We have used this version rather than 



The diagonalization of ( 12 ) for PBC in the limit of 



His — ~ z2 x * x i+i 



(12) 



due to the fact that having a diagonal interaction term 
is more convenient for the numerics. Of course both 
versions are equivalent since they can be transformed 
into each other by means of the unitary transformation 
U = i=1 Hi where the Hi are 1-qubit Hadamard gates. 
For the exact diagonalization however we have used ( 12 ) 



in order to stick closer to the procedure given in [!|| where 
the authors treat the spin-1/2 XY chain. 



an infinite number of sites is straightforward [10 . The 
first thing one has to do is to map the spin Hamiltonian 
to a fermionic one via a Jordan- Wigner transformation. 
Now the Jordan- Wigner transformation is non-local due 
to the fact that it transforms local spin operators into 
fermionic ones that anticommute if they act on different 
sites. Luckily for almost all terms in the Hamiltonian 
the non-localities cancel except for the term representing 
the interaction between the last and the first site. This 
term ends up containing a global parity operator acting 
on all sites and thus breaking the translational invari- 
ance with respect to the fermionic modes. Now, if we 
are interested in the thermodynamic limit, we will even- 
tually take the limes N — > oo at some point, and in this 
limit the contribution of one interaction term to the en- 
ergy can be neglected. We thus have the freedom to alter 
this term as we please in order to simplify things. One 
very convenient choice are the so-called Jordan- Wigner 
boundary conditions which are nothing more than sim- 
ply neglecting the global parity operator in the last term 
thereby yielding the fermionic Hamiltonian translation- 
ally invariant. Note that the Jordan- Wigner boundary 
conditions cannot be expressed in a trivial way in terms 
of spin operators. 

The fermionic Hamiltonian obtained in this way is 
quadratic and translationally invariant, but it is not par- 
ticle conserving. This can be fixed by a canonical trans- 
formation jS] to non-interacting Bogoliubov fermions. 
The ground state of the system is then given by the new 
fermionic vacuum while excited states can be obtained 
by sequentially filling the fermionic modes. Ordering the 
eigenstates of the original spin model by momentum and 
energy, it turns out that the lowest energy branch co- 
incides with the dispersion relation of the Bogoliubov 
fermions. This happens because for a given momentum, 
the lowest energy state is always a state where precisely 
one fermionic mode is occupied. 

Now, for finite systems with periodic boundary condi- 
tions, the Hamiltonian after the Jordan- Wigner transfor- 
mation presents a difficulty: due to the fact that in this 
case we cannot choose the boundary conditions freely, 
there is one term that contains a global parity operator 
as prefactor (see Eq. 2.11' in [9]). At first sight, this 
term makes the Bogoliubov transformation impossible. 
However, if we project the Hamiltonian onto the sub- 
spaces with either odd or even parity, we can replace the 
parity operator by its eigenvalue in that subspace s.t. it 
becomes ±1, and we can apply the Bogoliubov transfor- 
mation in each subspace separately. The spectrum of the 
original Hamiltonian can then be constructed by picking 
from each subspace the states with the correct parity. 
It turns out that we can arbitrarily choose the sign of 
the Bogoliubov parity by shifting the Fermi surface. For 
example, if we choose the fermionic vacuum to be the 
state with lowest energy, all excited states are particle 
excitations |11| and it turns out that the parity opera- 
tor changes its sign under the Bogoliubov transformation 
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FIG. 4. (Color online). Relative precision of the low excitation spectrum for a critical Ising chain with N — 50. Left: D = 8. 
Right: D = 32. 








1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


11 


12 


13 


14 


15 


16 


17 


18 


19 


20 


21 


22 


23 


24 


25 


1 


l») 


11) 


2 2' 


2 2' 


2' 2' 










l^,i> 
2 ' 2' 


|i?,i> 
2 2' 




|12) 


|13) 


|14) 


|16) 


|16) 


117) 


|18) 


|19) 


|20) 


|21) 


|22) 


(23) 


(24) 


125) 


2 


\-,-k 
2 2 


l 3 -,--> 


|2> 


13) 


l|-l> 


I-, 3 -) 

2' 2 ; 


\- 3 -) 

'2' 2' 




18) 


M 


HO) 


1") 


23 1 

1 2 2' 


25 1 
2 2' 




2 2' 


1 2 2' 


33 1 

33 1 

1 2 2 


\ 2 -.- 2 ) 


I 37 ,-) 
1 2 '2' 


39 1, 
'2*2' 


2 ' 2' 


I 43 ,-) 
1 2 ' 2 ; 




1^ 




3 


1,-1,0) 


2,-1,0) 


I 5 ---) 
2 2' 


12,1,0) 


|4) 


15) 


*!> 


1?) 


"¥■§> 


| 15 ,?> 
2 2 


I",-) 
2 ' l' 


I 18 , 3 -) 

2 2' 


i 21 .?) 

2 ' 2' 


1 12, 1,0) 


13,1,0) 


14,1,0) 


(15, 1,0) 


(16, 1,0) 


(17, 1,0) 


39 1, 
2 ' 2' 


| 41 , i) 
2 2' 


I 43 ,- 1 ) 
1 2 ' 2' 


45 1. 

2 " 2' 


| 47 ,- ! ) 
2 2' 


49 1 

2 ' 2' 


,-49 -1 
2 ' 2 ' 


4 


2' 2' 


if. §> 


3,-1,0) 


l-b 


13,1,0) 


14,1,0) 


w 




>".§> 


18,1,0) 


ft 1,0) 


|10, 1,0) 


111,1,0) 


23 3 

'T'2' 


2 ' 2' 


2 ' 2' 


1 2 ' l' 


1 2 ' 2 ; 


37 1 ,. 
1 2 


|18,1,0) 


119,1,0) 


|20,1,0) 


|21, 1, 0) 


|22, 1,0) 


|23,1,0) 


|24, 1,0) 



TABLE I. (Color online). Quasi-particle structure of the lowest four branches for g = 1. The red/blue (grayscale: dark/light) 
boxes highlight states from the odd-parity subspace respectively from the even parity subspace. The ground state which is not 
shown in the table is the fermionic vacuum in the even parity subspace i.e. \GS) = |^) e „ en - 
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FIG. 5. (Color online). Exact solution for the dispersion 
relation of the Bogoliubov modes at criticality (i.e., g = 1). 



for fields below the critical point i.e. g < 1. For g > 1 
this choice of the vacuum state leaves the parity opera- 
tor invariant. Thus for g < 1, in principle we can switch 
the sign of the parity operator by shifting the Fermi sur- 
face and thereby we could always define the Bogoliubov 
modes such that the parity operator remains invariant. 
We will however give numerical evidence for the fact that 
choosing the Fermi surface to be the fermionic vacuum 



state is the physical choice. 

Let us first present the results obtained at the critical 
point g = 1. In Figure [3] we have plotted the energy 
of the lowest ten branches of excitations of a chain with 
50 spins obtained for MPS bond dimensions D = 8 and 
D = 32. The results for D = 32 are so close to the exact 
spectrum that it makes much more sense to look at plots 
of the relative energy precision rather than at plots of the 
energy itself. This is shown in Figure |4j 

At first sight the crossing of the precision line for the 
first branch of excitations with the one for the second 
branch seems a little unusual. How can it be that states 
with higher energy are approximated by roughly two or- 
ders of magnitude better than states with lower energy? 
The answer to this question is obvious if one looks at 
how the eigenstates emerge from the elementary Bogoli- 
ubov modes. Tabic |l] shows which Bogoliubov modes 
contribute to each individual eigenstate in the first four 
branches of excitations. Modes from the even parity sub- 
space have half-integer momentum while the ones from 
the odd parity subspace have integer momentum. Note 
that since only excitations with an even number of par- 
ticles are allowed in the even-parity subspace, the result- 
ing states always have integer momentum. Henceforth 
\ty even shall denote the vacuum in the even parity sub- 
space and \Q) odd the vacuum in the odd parity one. The 
ground state of the critical chain is the Bogoliubov vac- 
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FIG. 6. (Color online). Relative precision of the low excitation spectrum for the critical Ising chain with different chain 
lengths. Left: TV = fOO, D = 32. Right: TV = 500, D = 20. 



uum in the even parity subspace i.e. \GS) = |^} eKen - 
The first excited state is the zero momentum state from 
the first branch and is given by a Bogoliubov mode with 
zero momentum from the odd parity subspace [12 . It is 
sufficient to show in Table [I] how the spectrum emerges 
from elementary excitations for momenta < k < TV/2 
since the dispersion relation of the Bogoliubov fermions is 
symmetric around k = as can be seen in Figure [5j The 
important thing to notice in Table [I] is that the lowest 
branch of excitations does not contain solely one-particle 
excitations as it does in the case of the infinite chain. 
Looking back at the right plot in Figure [4] we see im- 
mediately that the one-particle excitations from the first 
two branches are approximated with roughly the same 
accuracy between 10~ n and 10~ 9 with the lower value 
for small momentum states. One can easily check that 
the states with the same order of accuracy from higher 
branches are precisely one-particle excitations. On the 
other hand it is obvious that two-particle excitations from 
any branch where one of the particles has fixed momen- 
tum k = 1/2 can be found in the plateau with relative 
precision of roughly 10 -7 . The other plateaus of simi- 
lar precision in the D = 32 plot of Figure [4] represent 
either two-particle states where each particle has higher 
momentum or three and more particle excitations. 

This interpretation of the branch crossings in Figure [4] 
is strong evidence for the fact that Q is a very good 
ansatz for one-particle excitations. However it turns out 
that if D is large enough, ([!]) is also a fairly good ansatz 
for many-particle excitations. The reason for this is that 
the large bond dimension compensates for the localiza- 
tion of excitations inherent in ansatz ([!]) by spreading 
the effect of the optimized tensor B to a region around 
it whose size is of the order of the induced correlation 
length of the MPS we start with. This is exactly why 
for the Ising chain with g = 1, TV = 50 and D = 32 
even states from tenth branch are approximated with an 
accuracy of roughly 1CT 4 . 

Now let us have a closer look at the region of the 



"level-crossing" between the lowest two-fermion branch 
from the even parity subspace and the lowest one-fermion 
branch from the odd parity subspace. In the case of g = 1 
this crossing turns out to be at approximately TV/4. In 
the immediate neighborhood of the crossing the energy 
difference between states with identical momentum be- 
comes very small. Now if the bond dimension D is cho- 
sen such that the precision of the MPS is of the same or- 
der like the interlevel spacing, these two levels cannot be 
discriminated properly by the MPS algorithm and thus 
there is no clear interpretation we can give to these MPS 
states in terms of one or two-particle states. As can be 
seen in the D = 8 plot of Figure |4j in this region the first 
and second MPS branch interpolate between the one and 
the two-particle states which we can safely discriminate. 
Note that this observation holds only on the side of the 
level crossing where the one-particle state has higher en- 
ergy than the two-particle state (e.g. at k « TV/4 on the 
left side of the crossing). On the other side, the one- 
particle excitation has lower energy and our one-particle 
MPS ansatz is perfectly suited to discriminate between 
the first and the second branch even if the precision is 
smaller than the actual gap between the levels. 

The last thing we would like point out about Figure [4] 
is the gap in accuracy between the states from the second 
and third branch at momentum k = TV/2. It turns out 
that this is a doubly degenerate state since it can be cre- 
ated by two different superpositions of elementary excita- 
tions with the same energy namely | , \ ) and | — -y , — |) . 
This is the reason why the precision of the k = TV/2 
state in the second branch is better than that of the sur- 
rounding two-particle states which are not degenerated: 
variational algorithms are more precise if they try to ap- 
proximate the energy of an entire subspace of the Hilbert 
space rather than that of a single state. However since 
all states generated by our algorithm are orthogonal, the 
price we have to pay for the improved precision in the sec- 
ond branch, is a slightly worse precision of the k — TV/2 
state in the third branch. 
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FIG. 7. (Color online). Relative precision of the low excitation spectrum for the Ising chain at g = 1.1 for different chain 
lengths. Left: TV = 50, D = 32. Right: TV = 100, D = 32. 
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TABLE II. (Color online). Quasi-particle structure of the lowest three branches for g — 1.1. The red/blue (grayscale: 
dark/light) boxes highlight states from the odd-parity subspace respectively from the even parity subspace. The ground state 
which is not shown in the table is the fermionic vacuum in the even parity subspace i.e. \GS) = \Q,) even . 



With this said, we can present the results we have ob- 
tained for different chain lengths N and different values of 
the magnetic field g. Figure [6] shows the accuracy of the 
algorithm for chains with 100 and 500 sites at g = 1. The 
plot for TV = 100 is very similar to the D = 32 plot from 
Figure |4j At small momenta 6 < k < 11 the one-particle 
excitations lie in the branches 4 to 6. These states are 
not reliably reproduced by our algorithm within the pre- 
cision that is otherwise reached for one-particle excita- 
tions. Presumably this would be fixed by increasing the 
bond dimension D beyond 32. However at the moment 
we cannot go to larger D for N = 100. For N — 500 the 
maximally accessible bond dimension is D = 20. The 
corresponding plot from Figure [6] is very similar to the 
small D plot for N = 50. Again in the region of the 
" level-crossing" between one-particle and two-particle ex- 
citations around k — N/4 our algorithm has difficulties 
obtaining the maximally reachable precision for the one- 
particle states. 

Now let us look at how the algorithm performs when 
we move away from the critical point. Figure [7] shows 
the relative energy difference of the MPS approximation 
for g — 1.1 i.e. above the critical point. The most strik- 
ing feature in this regime is clear separation of the lowest 
branch of excitations from the higher ones. This happens 
due to the fact that in this case the lowest branch con- 
tains only one-particle states as can be seen in Table [Tl| 
Again if D is large enough (e.g. D = 32 for N = 50), 



the different plateaus of similar precision become clearly 
visible. The first one at roughly A re iEi(k) ~ 10~ 8 con- 
tains two-particle excitations from the second and third 
branch where one of the fermionic modes has momentum 
k = 1/2. The second one with precision around 10~ 6 
consists of states where one of the fermionic modes has 
momentum k = 3/2. Note that in the plot for N = 100 
the lowest branch has slightly better precision than the 
one in the N = 50 plot even though the virtual bond di- 
mension is the same. This happens presumably because 
in this case the chain is long enough such that the run- 
ning particle cannot "feel its own tail" due to the PBC. 
This is another piece of evidence that ansatz ([!]) is very 
well suited to describe one-particle excitations. Whether 
many-particle excitations are faithfully reproduced de- 
pends strongly on the magnitude of D with respect to 
N. 

For g < 1 the picture changes dramatically. We can 
see in figure|8]that at g — 0.9 the best precision for states 
from the lowest branch is five to seven orders of magni- 
tude worse than for g = 1.1. Without any knowledge 
of the quasi-particle structure of the spectrum this huge 
difference might look a bit surprising. Even more sur- 
prising is the fact that the best precision at g — 0.9 is 
one order of magnitude worse than at the critical point 
g = l. However looking at the quasi-particle structure 
in Table [TTT| clarifies the situation. As already mentioned 
above the parity of the Bogoliubov fermions in the odd- 
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FIG. 8. (Color online). Relative precision of the low excitation spectrum for the Ising chain at g 
lengths. Left: TV = 50, D = 32. Right: N = 100, D = 32. 
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TABLE III. (Color online). Quasi-particle structure of the lowest three branches for g = 0.9. The red/blue (grayscale: 
dark/light) boxes highlight states from the odd-parity subspace respectively from the even parity subspace. \Q) odd denotes the 
fermionic vacuum in the odd-parity subspace. The ground state which is not shown in the table is the fermionic vacuum in the 
even parity subspace i.e. \GS) = \Q) even - 



parity subspace can in principle be arbitrarily chosen by 
shifting the Fermi surface. Throughout this work have 
made the most natural choice of choosing all modes to 
have positive energy i.e. none of the quasi-particle exci- 
tations are hole modes. For g < 1 this choice switches 
the sign of Bogoliubov parity operator such that we must 
pick states with an even number of excitations from the 
odd-parity subspace. One might argue against this con- 
vention and claim that it would be much more natural to 
pick the Fermi surface s.t. the zero-momentum mode is 
a hole excitation which yields the Bogoliubov parity op- 
erator identical to the spin parity operator. In this case 
we would have to construct all states from this subspace 
using an odd number of quasi-particles. On the other 
hand Table |III| c learly shows that our one-particle exci- 
tation ansatzpj) is a poor approximation to all states in 
this regime thereby indicating that indeed for g < 1 there 
exist no one-particle excitations. Thus our choice of the 
Fermi surface is justified and we have to construct the 
spectrum by picking the even quasi-particle excitations 
from the odd-parity subspace. 

We can understand this behavior from another point 
of view if we consider an infinite chain with open bound- 
ary conditions. It is well known that in the region of the 
phase diagram where the ground state is doubly degener- 
ated, the elementary excitations are kink excitations. If 
we would however impose periodic boundary conditions 



on the infinite chain, the single kink states would not be 
eigenstates any more since the existence of one domain 
wall would automatically imply the existence of a sec- 
ond one. In finite systems with PBC, the situation is a 
bit more complicated since the ground state degeneracy 
is not exact (the energy difference decays exponentially 
with N), but we can still argue along similar lines that lo- 
calized perturbations, that interpolate between the states 
of the almost degenerated ground state manifold, must 
always come in pairs. 

B. Heisenberg model 

The other model we have studied is the antiferromag- 
netic (AF) Heisenberg spin-1/2 chain. The Hamiltonian 
reads 

JV N 
i=l i=l 

(13) 

where S a — a a /2 and a a denote as usually the Pauli 
operators. As we already mentioned, the tensors A that 
constitute the backbone of ansatz ([T]) are the results of 
the simulations presented in (5]. In that work we have 
obtained a TI MPS approximation of ground states for 
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finite spin chains with PBC using matrices Ai that were 
real and symmetric. These results themselves were based 
on previous work (T3] where we have approximated the 
ground state of infinite OBC chains by TI MPS with real 
symmetric matrices. Thus the starting point in the en- 
tire procedure that leads ultimately to the excited states 
presented here is the imaginary time evolution for an in- 
finite chain with a set of real symmetric matrix product 
operators (MPO). As we explained in [13J it is not possi- 
ble to construct these directly from the the Hamiltonian 
(13). However, by means of the unitary transformation 

U = [/^ — Yl^li crfj-i (i- e - ac ting with a o^-gate on 
every second site) we obtain 



Hi 



HB 



U^H HB U 



1 N 

i=l 



(14) 

which allows us to express the imaginary time evolution 
in terms of real symmetric MPO. Note that in order for 
this procedure to work we have to restrict ourselves to 
chains with an even number of sites. In this case it does 
not matter if we apply the a^-gates on sites with an odd 
or an even index, so without loss of generality we will 
apply them on the odd ones. Now Hhb and H' HB have 
the same spectrum and since their eigenstates are simply 
related to eachother by 



N/2 



2j 



-1 Wi) 



(15) 



we can digonalize H' HB first and obtain the eigenstates 
of Hhb subsequently with very little effort. 



X 





1 


2 


3 


4 


5 


6 


7 


8 


1 


— 8 


+1 


+2 


+3 


— 4 


— 3 


— 2 


— 1 


+8 


2 


— 8 


— 7 


+ 2 


+3 


+4 


— 3 


— 2 


— 1 


+8 


3 


+0 


— 7 


— 6 


-6 


—4 


+5 


+6 


+7 


—o 


4 


+0 


+ 1 


— 6 


~ 5 


+4 


— 3 


+6 


+7 


—o 


5 


+0 


+ 1 


— 6 


+ 3 


+4 


— 3 


+6 


— 1 


—o 


6 


+0 


+ 1 


— 6 


—5 


—4 


+5 


+6 


— 1 


—o 


7 


— 8 


1 1 


— 6 


—5 


+4 


+ 5 


+ G 


— 1 


+8 


8 


— 8 


— 7 


— 6 




—4 


+ 5 


+C 


— 1 


+8 


9 


+0 


— 7 


+ 2 




+4 


+5 




+7 


—o 


10 


+0 


— 7 


+2 


+3 


—4 


+5 


— 2 


+7 





TABLE V. (Color online). Multiplet structure of the low- 
est ten branches of excitations for a Heisenberg 16-site chain 
with Hamiltonian \14) . The colors encode the multiplet infor- 
mation: yellow-singlet, blue-triplet, red-quintuplet (grayscale: 
darker colors encode higher multiplets). The sign denotes the 
parity of a state and the index denotes the momentum k if 
we apply the transformation ( 15 1 to an eigenstate with mo- 
mentum k! . 



We will show below that the momentum of a state is 



not always invariant under the transformation (14). The 
easiest way to obtain the momentum for any given state 
is by computing the expectation value of the translation 
operator T with respect to that state. Hhb and H' HB are 
both translationally invariant thus all their eigenstates 
have well defined momentum so we can be sure that the 
reverse transformation |?//(fc')) — > \ipi{k)) will map mo- 
mentum eigenstates to momentum eigenstates albeit k 
will generally differ from k' . The relation between the 
momenta follows easily from 
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TABLE IV. (Color online). Multiplet structure of the low- 
est ten branches of exc itations for a Heisenberg 16-site chain 
with Hamiltonian (13 1. The colors encode the multiplet in- 
formation: yellow-singlet, blue-triplet, red-quintuplet, dark 
red-septuplet (grayscale: darker colors encode higher multi- 
plets). The states within each multiplet are ordered according 
to their total spin projection quantum number. The sign de- 
notes the parity of a state. 



e -<-S" = (^(k)\T\^(k)) = 

N/2 N/2 

= (#')! 114-iW 114-1 MW) 

S=i ' S=i ' 

N/2 N/2 

= mk')\ IK--i IK mm')) 



3=1 



N 



= e"^ mk')\ n o* mk')) - (P y ) il>v 

3=1 

(16) 

where we have used T^O^T- 1 = l\j°3+i and 

T|^(fc')) = Wi{k')). Thus the change in mo- 

mentum depends solely on the expectation value of the 
operator P y — YljLi a% j which in the following we will 
call the 'parity operator. This naming convention makes 
sense since P y — i N exp^nS^) where = Ylj=i $j 
thus P y measures the parity of the total spin along the 
y-direction. Note that due to the factor i N the mean- 
ing of positive and negative parity is interchanged for 
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chains with N — 0( mod 4) and chains with N = 2( 
mod 4). The parity is a good quantum number for 
both H HB and H' HB so there exist eigenstates ^[(k 1 )) 
that have well defined parity plus or minus one. If 
= +1 the momentum remains unchanged i.e. 

we have k = k' ®n N/2 



P 



V/i'k' 



k = k', if (Py 



i'k' 



= -1 = e 



where (Bn denotes addition modulo N, Note that the 
parity itself is invariant under the mapping between Hjj B 
and H' HB since W P y U — P y . 

Now the generators of the SU (2) symmetry for H' HB 
do not commute with the translation operator thus we 
cannot classify the momentum eigenstates in terms of ir- 
reducible representations of SU(2). For Hhb however 
we can do this, so we know exactly the degeneracy struc- 
ture of the spectrum in each subspace with fixed mo- 
mentum. Thus if we encounter for instance a three- 
fold degenerated eigenstate of H' HB , we know this is 
mapped to a spin triplet with well defined momentum 
in the original Hamiltonian. Accordingly it must contain 
a two-dimensional subspace with negative parity corre- 
sponding to total spin along the y-direction ±1 and a 
one-dimensional subspace corresponding to total spin 0. 
Since the spin triplet in the original Hamiltonian has well 
defined momentum, according to the rules for the map- 
ping k k', we will have one eigenstate of H' HB with 
momentum k and a two-dimensional subspace with the 
same energy but different momentum k' = k Qpf N/2. 
In this way |15j . after approximating the spectrum of 
H' HB and labeling all energies with the corresponding 
momentum we can obtain the spectrum of Hhb by mere 
inspection of the degeneracy structure. Table [TV] and Ta- 
ble [V] illustrate how the multiplets of Hhb and H' HB are 
related to eachother. 

This procedure works very well for the lower branches 
of the dispersion relation where the precision of our sim- 
ulation is good enough to discriminate unambigously be- 
tween different multiplets. For higher branches, on one 
hand the precision gets worse and on the other hand the 
density of states increases such that multiplets with sim- 
ilar energy become effectively undistiguishable for our 
algorithm. In this case the eigenstates with well defined 
momentum that we obtain for the Hamiltonian H' HB do 
not have well defined parity i.e. they mix parity eigen- 
states with different parity. Since states with same mo- 



mentum and different parity are mapped by ( 15 ) to states 



with different momentum, if we start with such a mo- 
mentum eigenstate we obtain after the transformation a 
superposition of states with different momenta which is 
clearly not a momentum eigenstate. There are however 
two ways to overcome this issue and obtain approxima- 
tions of the eigenstates of Hhb that are at the same time 
exact momentum eigenstates. 

The first one amounts to computing the matrix el- 
ements of the translation operator T in the subspace 
spanned by the transformed states {M^ dd \il>i(k)')} where 

M^ dd := rifci^Ij-i ar± d then diagonalize this ma- 
trix. It is not difficult to check that this can be done 
for each momentum k' separately since M v odd T M^ dd = 



^odd Meven T = P y T which is a translationally invari- 
ant operator and thus it does not mix states with dif- 
ferent momentum. Diagonalizing each of the T(k / ) i j = 
U> (k')uD(k')i m V (k') m j yields for each k! a unitary U (k ! ) 
that is nothing more than the transformation that we 
need to obtain the desired momentum eigenstates via 
IV^(^i)) = U{k')ij fy'Ak 1 )). The new momentum ki can 
be read off the diagonal matrix D(k'). There are two 
drawbacks that come with this procedure. The first one 
is that we must compute the matrix elements T(k')ij each 
of which is done with the computational cost 0(ND 5 ). 
Since there are Nb 2 of these where b is the number of 
branches, we obtain the overall cost 0(N 2 b 2 D 5 ). Usu- 
ally we compute enough branches such that b 2 > D holds, 
thus the cost for this procedure ends up being higher 
than the one for the diagonalization of H' HB . The second 
drawback is that the superpositions U(k')ij \ip'j(k')) mix 
the original approximations of the energy levels thereby 
slightly lowering the energy of higher excitations but in- 
creasing the energy of lower excitations, which are usually 
the ones we are most interested in. 

The second way to approximate the eigenstates of 
the original Hamiltonian Hhb such that they are at 
the same time exact momentum eigenstates is to add 
to H' HB a perturbation that splits degenerated levels 
with different parity. This is easily achieved by taking 
H BB := H HB ±XP y where A must be chosen such that it 
is big enough for our algorithm to deliver only states with 
a single parity, but as small as possible in order to avoid 
numerical inaccuracies caused by altering the Hamilto- 
nian. In the case of the Heisenberg model, if we choose 
to compute b = 10 branches, A = 0.1 • N fulfills these 
requirements. In practice we first apply our algorithm to 
H HB which yields for each momentum k' b states with 
positive parity. These states do not change their mo- 
mentum under the transformation ( |15[ ). Subsequently 
we apply the algorithm to H BB which yields states with 
negative parity that change their momentum after the 
transformation according to k = k' ®n N/2. In this way 
we end up with 2b branches of states that approximate 
the spectrum of Hhb and that are at the same time 
exact momentum eigenstates. The computational cost 
per state is thus exactly the same like diagonalizing only 

h 'hb- 

Let us first look at the results we have obtained for 
a small chain with 16 sites. We have chosen to look 
at such a small system first for two reasons: First, 
even though the Heisenberg model is exactly solvable via 
Bethe ansatz, obtaining all energy levels can be quite in- 
volved. Choosing N as small as 16 allows us to compute 
the spectrum of this small chain by means of exact diag- 
onalization. Second, even for the energy levels that are 
easily computable with the Bethe ansatz (i.e., the triplet 
states in the subspace of two-spinon excitations [5]) it is 
not possible to obtain the eigenstates themselves. Exact 
diagonalization of a small chain on the other hand allows 
us to compute and store the exact eigenstates in order to 
check the fidelity of our MPS approximation. 
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FIG. 9. (Color online). Results for the low excitation spectrum (left) and the corresponding relative precision (right) for the 
Heisenberg spin-1/2 chain with TV = 16 sites. 
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FIG. 10. (Color online). Distance between several degener- 
ated subspaces obtained by our algorithm to the correspond- 
ing degenerated subspaces obtained by exact diagonalization. 
As a measure for the distance we have used the sine of the 
canonical angle with the largest magnitude as defined in [14] . 



Figure [9] shows the energy of the first ten branches 
of excitations and the corresponding relative precision. 
Note how states belonging to the same multiplet have 
very similar precision even though they have different 
parity and thus correspond to eigenstates of H' HB with 
different momentum. Since there are no one-particle ex- 
citations in the low-energy spectrum of the AF Heisen- 
berg model, we do not obtain such a good precision like 
in the case of the quantum Ising model. Nevertheless we 
get a very good approximation of the first excited level, 
namely the triplet excitation at k = N/2. We have also 
tested the accuracy of the states themselves: for non- 
degenerated states, the absolute value of the overlap is 
a perfect measure for this, and for reasons that will be- 
come clear immediately, we have looked at the sine of the 
fidelity. For degenerated states, in order to compare the 
subspace spanned by our MPS to the one spanned by the 
exact eigenstates, we have used as a measure for the dis- 



tance the definition given in chapter 7 of [2] : the sine of 
the largest canonical angle between the two subspaces. 
The canonical angles can be easily computed from the 
matrix that has as its entries the overlaps between all 
states of the subspaces that we want to compare. The 
results are plotted in figure[l0] We see that only the MPS 
with momentum k = and k — N/2 are extremely accu- 
rate. All other states, especially those around k = N/4, 
are much further away from the exact solutions, which 
is a bit surprising given the fact that the energy preci- 
sion for these states is comparable to the one obtained 
for k — 0. 

The spectrum that we obtain for longer chains is plot- 
ted in figure [TTJ In this regime we have only looked at 
the precision of the lowest two-spinon triplet for which 
the exact results were obtained following [8]. Again we 
see that the states at momentum k = k ® N/2 have the 
best accuracy. We would like to make two more remarks 
concerning the chain with N = 50. First, note that the 
ground state has momentum ko = N/2 in this case. Sec- 
ond, unlike in the case of N = 100, where for all momenta 
k ko the lowest excitation is a triplet, we observe that 
for N = 50 this does not happen. Our simulations reveal 
that at k G {2,3,47,48} the quintuplet excitation lies 
below the triplet while at k € {23, 27} it is a singlet that 
is the lowest lying excitation. 



V. CONCLUSIONS AND OUTLOOK 

Inspired by previous approaches [TJ |3] we have intro- 
duced a new method for the simulation of translationally 
invariant spin chains with periodic boundary conditions. 
We have used an MPS based ansatz that corresponds 
to a particle-like excitation with well defined momentum 
in order to obtain extremely accurate results for models 
where the spectrum contains precisely one-particle states. 
For states that can be expressed in terms of many quasi- 
particle excitations, we still obtain feasible results if the 
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FIG. 11. (Color online). Results for the low excitation spectrum and the corresponding relative precision of the lowest triplet 
(insets) for the Heisenberg spin-1/2 chain with N = 50 (left) and TV = 100 (right) sites. 



MPS bond dimension is chosen to be big enough. In the 
case of the quantum Ising model, our results indicate that 
for g < 1 the spectrum is built up entirely out of exci- 
tations with an even number of quasi-particles. General- 
izations of our approach can go in two directions: First, 
it is possible to adjust ansatz |l]) in order to treat infinite 
systems with open boundary condition, which we are ad- 
dressing in 7\ . Second, it seems feasible to generalize our 
approach to a many-particle ansatz by using more than 
one MPS tensor in (JlJ in order to define the variational 
manifold. 
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